function ovecs = orthvecs(ivecs, vorder, detr)
% orthvecs  - orthogonalize vectors against each other
%
% FORMAT:       ovecs = orthvecs(ivecs [, vorder, detr])
%
% Input fields:
%
%       ivecs       NxV matrix with V number of vectors
%       vorder      must contain all numbers from 1..V, default 1:V
%       detr        if given and evaluates to true, detrend vectors
%
% Output fields:
%
%       ovecs       orthogonalized vectors
%
% Note: for two vectors, the "cheaper" orthvec function can be used

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ~isa(ivecs, 'double') || ...
    isempty(ivecs) || ...
    ndims(ivecs) > 2 || ...
    size(ivecs, 1) < size(ivecs, 2) || ...
    any(isinf(ivecs(:)) | isnan(ivecs(:)))
    error( ...
        'BVQXtools:BadArgument', ...
        'Valid matrix of vectors without Inf/Nans are required.' ...
    );
end
numvals = size(ivecs, 1);
numvecs = size(ivecs, 2);
if nargin > 1 && ...
    isa(vorder, 'double') && ...
    numel(vorder) == numvecs && ...
    numel(vorder) == max(size(vorder)) && ...
    all(vorder > 0) && ...
    numel(unique(fix(vorder))) == numvecs
    vorder = fix(vorder(:)');
else
    vorder = 1:numvecs;
end
if nargin > 2 && ...
   ~isempty(detr) && ...
   (isnumeric(detr) || islogical(detr)) && ...
    detr(1)
    detr = true;
else
    detr = false;
end

% generate output
ovecs = zsz(ivecs);
tcol = vorder(1);
if detr
    ovecs(:, tcol) = ivecs(:, tcol) - mean(ivecs(:, tcol));
else
    ovecs(:, tcol) = ivecs(:, tcol);
end

% iterate in order of importance
for vc = vorder(2:numvecs)
    
    % get vector which is to be tested
    testvec = ivecs(:, vc);
    if detr
        testvec = detrend(testvec);
    end
    
    % generate regression matrix
    rmat = ovecs(:, tcol);
    rmat = rmat - repmat(mean(rmat), [numvals, 1]);
    tmat = rmat';
    frmr = pinv(tmat * rmat) * tmat;
    
    % regress
    b = frmr * testvec;
    
    % fitted -> residual
    r = testvec - rmat * b;
    
    % residual goes into output !
    ovecs(:, vc) = r;
    tcol = [tcol, vc];    
end
